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Abstract - The goal of this paper is to set up a framework designed to take into account the character- 
istics of sediment particles when transported by water. Our protocol consists in describing the characteristics 
of sediment particles via an additional variable, and to build operators involving this new variable, modeling 
the evolution of the particle characteristics. Several such operators are proposed, some based on principles of 
relaxation toward an equilibrium, and others on a description of the particles' aggregation and fragmentation 
process. A discrete version of the latter is also offered for numerical settings. 
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1 Introduction 



In view of the evolution of the climate, and the increasingly stringent requirements in terms of 
feasibility and impact studies before dredging or building sea walls or harbors, the estuary morpho- 
dynamical issue is becoming a topic of major significance. As a consequence, the behavior of the 
complex sediments found in estuaries is today widely studied, analyzed and modeled. The research 
effort concerns every aspect, from measurement protocols to simulations of deposition, erosion, 
transport by water, wave action, turbulence results and flocculation processes. 

In all these aspects, modeling has a key role to play. To summarize, the modeling of cohesive 
sediments involves three compartments, as well as the interactions which link them together. The 
first compartment deals with fluid field forecasts, and involves Navier-Stokes or Shallow Water 
type equations, possibly involving turbulence and using a propagating eddy viscosity. The second 
describes sediment behavior when it is deposited on the seabed, and the third models the transport of 
sediment particles when suspended in the water column. Roughly speaking, this last model propels 
sediment particles at the same velocity as the fluid added to a settling velocity. It may also take 
into account the action of turbulence, using the eddy viscosity as a dissipative effect on sediment 
particles. 

The present paper is situated in this context, and focuses on the issue of the transport of sediment 
particles by the water column. It offers a robust and flexible framework which takes into account the 
evolution and alterations of the sediment particle characteristics (size, mass, porosity, etc.) while 
the particles are in the water column, concomitantly with other phenomena (transport, settling, 
turbulence). The main idea consists in introducing a mass density of sediment particles p{t,'x.,X) 
depending on time t, position x = {x,y,z), and also an additional variable A € A, which describes 
the particles' characteristics. 

To present this idea, this article will begin by giving a proper definition of p, in the context of 
this framework, consisting in a two-part integro-diffcrcntial equation. The first part consists in a 
differential operator acting on p and describing the action of water as it transports particle. The 
second part is an integral operator, named G, modeling the evolution of the particles' characteristics. 
Examples of sets A of characteristics, variables A and integral operators G will then be given, 
before exploring the way in which existing aggregation models may be translated into the present 
framework. Finally, a discrete instantiation of set A and operator G are given for numerics. 

2 Guiding ideas 

The following elements will be taken for granted: A given estuary may be represented by a regular 
subset il e K'^, provided with coordinates x,y,z, where the x— axis is horizontal and points toward 
the east, the y— axis is horizontal and points toward the north, and the z— axis is vertical and points 
toward the sky. At any time t G and in any point x = (x, y, z) of il, the water velocity v with 
coordinates (ii,w,w) in (x, z)— coordinate system may be computed using a Navier-Stokes-type 
system, possibly involving eddy viscosity. The salinity S and the temperature T may be obtained 
by solving advection-diffusion equations, possibly involving eddy viscosity. Moreover, the energy of 
turbulence k and its dissipation rate e may be computed using, for instance, a k — e model. The 
eddy viscosity involved in the equation describing the evolution of v, S and T may be computed 
from k and e. We may also suppose that the water pH and the amount of organic matter per liter 
of water O are available. In the sequel, T denotes the fluid field 

^=(v,S,r,k,£,pH,0), (2.1) 

which ranges in x (M+)^ and which depends on t and x. 
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The main idea to be explored hereafter consists in assuming that the characteristics of the 
sediment particles may be described by a variable A belonging to a given continuous space A and 
that, at time t £ R"*", the mass distribution of suspended matter of type A G A and in point x G 
may be described by a measure, which is absolutely continuous with respect to the Lebesgue measure, 
and with density p(i,x. A). The precise definition of p states that for any subset uj C il x A, the 
mass of sediment particles with position and characteristics situated within lo is 



/ p{t,x,X)dxdX, (2.2) 

J UJ 



at any given time t. 

Actually, particles result from the assembly of elementary sediment particles. Hence, the char- 
acteristic of a given particle naturally belongs to a discrete space A, and the mass distribution R, 
with respect to the characteristic variable A, is naturally a sum of Dirac mass 

i? = ^A^(t,x,A)dx(5^^3^. (2.3) 
AeA 

Hence, when making the above assumption, we consider that, on the observation scale, R may be 
replaced (or approached) by p{t, x. A) dxdX. 

From p, the mass density of suspended matter at t is defined by 

r(t,x) = / p{t,x,\)d\. (2.4) 
Ja 

The evolution of p, over time, is supposed to be the result of transport by water, settling, diffusion 
by turbulence and aggregation, fragmentation and, more generally, the shape and mass evolution of 
the particles. In point of fact, p should be seen as the solution to: 

|..(..A,|.n^,X,|.(,n^,X,-«.,A,,)| 



a, ay ■ a. . = c:(^-aA)^ (2^5) 

The first four terms in (|2.5p are the time derivative of p following the trajectories induced by velocity 
([/, V,W- Ws). Velocity V = A) = (C/(JP, A), y(JP, A), WiJ", A)) is the velocity transmitted by 

water to the particles. It depends on the fluid field T and on the particles' characteristics A. But in 
most situations, it is reasonable to set V(J^, A) = v, meaning that water transmits its velocity directly 
to sediment particles. Settling velocity Ws(A, r) is the velocity at which particles fall toward the 
seabed. It is natural to consider that Ws depends on particle characteristics, with the idea that the 
heavier a particle, the faster it falls. Moreover, if the particle density r in the water column is high, 
settling may be slowed or hindered by the proximity of many particles. In this case, 144 has to depend 
on r. The fifth term of the equation's left-hand side conveys the fact that sediment particles undergo 
diffusion. This involves a horizontal diffusion coefficient /i and a vertical diffusion coefficient v. The 
diffusion phenomenon essentially comes from turbulence. Hence, choosing for /i(^. A) and v{J-, A) the 
water's usual eddy viscosity, given by c^— (where ~ 90) is not completely unreasonable, at least 
for particle characteristics A corresponding to small sizes. For particles characteristics corresponding 
to sizes bigger than the turbulent structures, other choices must be made. The right-hand side of 
eauation (|2.5p models the evolution of the particle characteristics. Clearly, the way in which particles 
combine, fragment, grow or, more generally evolve is linked with the fluid field, especially aspects 
such as temperature, salinity, pH, concentration of organic matter and certainly turbulence energy. 
Thus G depends on J^. The next section will provide examples of operator G and a discussion of 
its properties. 
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3 Relaxation models 



This section gives simple examples of operator G. 

In the following paragraphs, the operator G is not built up from physical considerations, but 
only by considering the asymptotic evolution of sediment particles when environmental conditions, 
or in other words, when the fluid field remains the same over a long period of time. 



3.1 A mass-preserving relaxation model with one-dimensional A 

In this first example, A is assumed to be and A G A stands for particle size. Considering this 
certainly supposes that all the sediment particles under consideration have the same shape, and that 
they can be characterized by a one-dimensional parameter. If the particles are onc-dimensional, A 
is the particle length, and particle mass is in direct proportion to A. If they are two-dimensional, A 
is the particle diameter and particle mass is in direct proportion to A^. If they are tridimensional, 
A is also the diameter but particle mass is in proportion to A'^. 



3.1.1 Operator building 

The operator-building process presented here is influenced by Bhatnagar, Gross & Krook 1 , who 
offered a kinetic model for gas dynamics. 

If it is well established that, when fluid field J- remains the same in a given place over a long pe- 
riod, the mass density distribution with respect to A is given by the following equilibrium distribution 
function: 

vT>eq(T,\), (3.1) 

where p is a mass density with respect to x— variable (p depends on x) and where 2?eg(.^i-) is a 
density probability defined by A (in particular, it satisfies: 

/ X'eq(^,A)dA = 1, foraU JP) (3.2) 

J K 

then, introducing a relaxation time 7^^, G may be defined as 

G(^,p,A) = -^(^p- (^pdA')Peg(^,A)^, (3.3) 



or 



G(.F(t,x),p(i,x,.),A) = -^(^p(i,x,A) - (^p(t,x,A')(iA')2?eg(^(i,x),A) 
3.1.2 Example of function Vf-q 

As function 'Deq{J-, A), defined for any G K."^ x (IR+)^ and any A G M^, we can choose: 
I?eg(.F,A) ^ ifA<A„,i„ 



(3.4) 



A — Amin /A — A,„iii\ . (3.5) 

^^P / /-7-\\ otherwise , 



which is drawn in figure [1] for Amin — 5 and for cr(J^) = 1 on the left and a{J-) = 3 on the right. 
This choice makes it possible to take into account that particles cannot be smaller than A,„i,i, and 
that the variability of particle size, for a given equilibrium, depends on a function cr(J-') of the fluid 
field. 
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Figure 1: Function Peq(-^, A) defined by p.Sp for Amin = 5, (t{T) — I (left) and a{!F) = 3 (right) 
3.1.3 Properties 

When used in (|2.5p . the operator G defined by (|3.3p pushes p{t,x,X) toward ( /9(t, x, A') dA') 



'Deqi^', A) at any time and place, with a relaxation time of Tlq. Moreover, the action of G in (|2.5 
does not influence the evolution of the total mass of sediment. 

In order to be more precise, a function p{t, A), not depending on x, which is solution to 



for a fixed vector G 



)'^, has the following properties. First, the quantity 



p{t, A) dX, 



remains constant over time and, secondly, for every A G A, the quantity 

p{t,X)-[j^p{t,X')dX')v,,{T, A), 
is divided by e after any period of time of length T^q. 
Property p.7p may be seen by integrating (|3.6p : 



(3.6) 



(3.7) 



(3.8) 



pdX 



dt 



A dt 



Teq J A 



p[.,X')dX')v,q{T, A)jdA 

^ p(.. A') dA - ^ p(.. A') dX^ = 0. (3.9) 



Property (|3.8p may be seen by computing the solution to equation (|3.6p leading, for any t > s, to 
p{t,X) - (^j^p{t,X')dX') V,q{J^,X) = [p{s,X) - (^p(s,A')dA') Peg(.F, A)^e(^-*)/^'. (3.10) 

3.2 A non-mass preserving relaxation model with one-dimensional A 

In cases when sediment cohesion is insured by a biological factor with an impact on the particles' 
mass, it is not reasonable to use a mass-preserving model. It is preferable to use a model able to 
reproduce the fact that for two sediment-particle populations issued from the same initial sediment- 
particle population - one made up of small particles and the other of large particles - the total 
mass of the second population is greater than the total mass of the first. Figure [5] shows two mass 
distributions with respect to A. Their total mass is not the same. A non-preserving mass model will 
be able to generate mass distributions of those shapes, from the same initial mass distribution. 
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Figure 2: Left: example of a small-particle population mass distribution. Right: example of a 
large-particle population mass distribution. 



Elementary sediment particle 
'^bio / Biological particle 



Figure 3: Example of particle resulting from the aggregation of 6 elementary sediment particles and 
5 biological particles. 



3.2.1 Operator building and properties 

To achieve the goal presented above, a decreasing function /(A) is introduced. Then, using 



1 



p(t,x, A')/(A')dA' 



Gi^it, x), X, .), A) = -i^[p(.t^ x> A) - ^ — n,iHt, x), A) ) (3.11) 



in (|2.5p pushes the solution toward 

p(i,x,A')/(A')dA' 



/(A) 



■Pe,(^(i,x),A), (3.12) 



with a relaxation time X,q and does not influence the evolution of 

[ p(t,x,A)/(A)dA. (3.13) 

3.2.2 Example of function / for one-dimensional particles 

A suitable function / for one-dimensional particles aggregated by means of biological factors may 
be built as follows. If the particles arc the result of assemblies of elementary sediment particles with 
length Amin, joined together with biological particles of length Abio, the length of a given particle is 

A = nAmin + (n- l)Abio, (3-14) 

for a given n. An example of such a particle, with n = 6, is represented in figure [31 Obviously, n 
can be expressed in terms of A, Amin and Abio- Indeed, since (rt — l)Aniin + (n — l)Abio = A — Amin, 
the following formula are true: 

A — Amin , A + Abio /oirA 

n — 1 — — and n = — . (3.15) 

Amin 4^ Abio Amin ^ Abio 
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The following quantities may also be computed: 

Q^^'-j _ '^-^min _ (-^ + Abio)-^min ^ ^ 

A (Amin + Abio)A ~ 

(1 — 0(A)) — (" ~ l)Abio _ (A - Amin)Abio ^ 
A (Amin + Abio)A 

Those quantities are the lineic proportions of elementary sediment particles and of biological particles 
in a particle of length A. As a matter of fact, the mass m{X) of a particle of length A may be expressed 
as 

m(A) = 0(A)AXmi„ + (1 - 0(A))AA^bio - ^^r^^^^^^-^-n + ^^r^^f^A^bio, (3.17) 

where A^min is the lineic mass density of elementary sediment particles, and A^bio is the lineic mass 
density of biological particles. 

Beside this, when a particle of length A joins up with another of length A', the result is a particle 
of length 



which satisfies: 



since 



A" = A + A' + Abio, (3.18) 
e{X")X" ^9{X)X + 9{X')X', (3.19) 

(A" + Abio)Amin (A + Abio)Amin (A' + Abio)Aniin 



(Amin + Abio) (Amin + Abio) (Amin + Abio) 

A/Iultiplying (|3.19|) by Afmin and rewriting it as 

0(A")A"A(min 

-m(A ) = 



(3.20) 



0{X")X"Mmin + (1 - e(A"))A"A1bi. 



^(A)AMnin 0(AOyAimin „^ . 



e(A)AAJmin + (1 - e(A))AMbio ' ' 0(A')A'Mmin + (1 - 0(A'))A'A^bio 

it may be deduced that defining 

, . ^ ^(A)A^n.in ^ 1 ^ 1 .o 02) 

•'^ ' 6'(A)Xmin + (1 - 6'(A))Aibio 1 - 6>(A) Atbio ^ A - Amin Abio-Mbio ' ^ • ' 

61(A) M mill 

A + Abio AminA^ mill 

equality 

/(A")m(A") = /(A)m(A) + /(A')m(A'), (3.23) 

is true. 

As a consequence, if a sediment-particle population is made, for each n G N, of A^i(A(ri)) particles 
of length A(n) = nAmin + {n — l)Abio, and if a second population is generated by joining and breaking 
down particles of the first population, leading to A^2(A(n)) particles of length A(n), for each n £ N, 
then the following link between sequences {Ni{X{n)y)n£n and (A^2(A(n)))„gN holds 

Ni{X{n)) f{X{n)) m{X{n)) = ^ A^2(A(n)) /(A(n)) m(A(n)). (3.24) 
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Since the mass measure, on A, of the two sediment-particle populations are 

Ri = Y.N, (A(n)) m(A(n)) , (3.25) 

r!6N 

R2^J2 ^2(A(n)) m(A(n)) 5a=a(„), (3.26) 

where 5\^x{^n) stands for the Dirac mass located in A(n), formula (|3.24p expressed also as 

<i?i,/>=<i?2,/> . (3.27) 

Now, if measures Ri and R2 are replaced (or approached) by pi(A) dX and P2(A) dX, formula (|3.27|) 
yields 

/ pi(A)/(A)dA= / p2(A)/(A)dA. (3.28) 

J A J A 

Then property p.l3p leads to the conclusion that choosing / defined by (|3.22p in the definition 
of G induces a behavior of p with respect to A in accordance with sediment that aggregates 
because of biological factors as described above. 

In cases when A^bio < -Mmin and Abio << Amin, since A > Amin, it is also clear that Abio << A. 
Then, /(A) may be approached by: 

A — Amin AbioA^bio ^ A — Amin AbioA^bio 29^^ 

A + Abio Amin A( min A AminAlmin 



3.2.3 Example of function / for any-dimensional particles 

For particles with dimension d (which may not be an integer) a suitable function may be built using 
a generalization of the above considerations. As previously, the considered particles are the result 
of assemblies of elementary sediment particles stuck together with biological particles. A particle 
with characteristic length A is considered to be made of a proportion 0(X) of elementary sediment 
particles, and of a proportion 1 — 6{X) of biological particles. The mass to(A) of such particles is 
then 

m(A) = 0(A)A'^A(min + (1 - e(A))A'^A(bio, (3.30) 

where Almin and Mmo are linked with a d— dimensional mass density of elementary sediment par- 
ticles and a d— dimensional mass density of biological particles. When a particle of characteristic 
length A joins with another of characteristic length A', it generates a particle of characteristic length 
A" such that 

^(A'OA'"^ = eiX)X'^ + d{X')X"^, (3.31) 

0(A")A"'*A(min 



0(A")A"'^Xmin + (1 - e{X"))X"''M^i 



-m(A") 



e(A)A'^A(min .y,^ 0(V)V'^.Mmin 



0{X)X<iMmin+{l-d{X))X'iMbio e{X')X"^Mmin+ {I - 0{X'))X"^Mb 

Then choosing 

f(X\^ ^(-^)-^min ,„„„^ 

^^''^ 0(A)A^min + (1 - e(A))A^bio ' ^ ^ ^ 
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insures the following 



/(A")m(A") = /(A)™(A) + /(A')m(A'), (3.34) 

leading to the conclusion that choosing /(A) defined by (|3.33p induces a behavior of p in accordance 
with the behavior of the distribution of d— dimensional particles, which aggregate by means of a 
biological factor. 



4 Statistical aggregation and fragmentation models 

In this section, aggregation and fragmentation models usually used in chemical engineering and 
colloid sciences and often referenced as "Monte-Carlo Simulation" (see Gardner & Theis [8], Spilman 
& Levenspiel , Daniels & Hughes [4| , Meakin [T^ , Liffman TT] , Shah et al. [l^ , Das [5] , Spouge 
[16], and Van Peborgh & Hounslow [iTj) are adapted to the framework presented in this paper. 



4.1 Statistical aggregation and fragmentation considerations 

Summarizing the ideas used in the above-cited references, and revisiting them with a viewpoint 
inspired from the Boltzmann equation context (see Cercignani [3]) we are led to the following rea- 
soning. The variable A € A = [Amin, +oo) stands for the particle size which is minimized by the size 
Amin of elementary sediment particles, and m(A) for the mass of particles of size A. If the particles 
are d— dimensional, 

m(A)=A/^dA^ (4.1) 

for a constant Md depending on d. 

Aggregation is described introducing a transition probability Ba{T , A, A'), which depends on the 
fluid field T. By definition, Ba{J-, A, A') is the probability that two particles, one of size A and one 
of size A', being in the same place in fluid conditions T , aggregate in a unit of time. It has of course 
the following property 

Ba{T,\,\')=Ba{T,\',\), (4.2) 

for any T ^ A and A'. (Expressions of transition probability Ba^ based on physical principles, may be 
found in Kim et al. |9j.) 

Once aggregated, the two particles give a particle of size A" = (A"^ + \"^Y/'^ which is such that 

m(A") = to(A) +m(A'). (4.3) 

Fragmentation is described using Bf{!F, A), which is the probability that a given particle of size A, 
in a place with fluid conditions J-', fragments in a unit of time. Then, Be{J-^, A, A') is the probability 
density function, with respect to variable A', that a particle of size A which fragments gives a particle 
of size A' < (AV2)^^'' (and another of size A" = (A'' - X"^)^/'' > (A''/2) ^^'^). By deflnition. Be has 
the following properties 

Be(.^,A,A') = Oif A' > (AV2)'^^ 
/ Be{T,X,X')dX' ^ [ Be{T,X,X)dX' ^l. ^^■^'^ 

Jx'eA J\'<[\''/2) 

Denoting by Be{J-, A, A") the probability density function, with respect to variable A", that a particle 
of size A which fragments gives a particle of size A" with (A''/2)^^'' < A" < A, and by L(A') = 



9 



(A'' - X"^Y^'^, L-\X") = {X'^ - A'"')!/'', for any set w C [(A''/2)^^'', A], Be and Be are linked by 

/ Be{J',X,X")dX" I Be{J',X,X')dX', (4.5) 



since every time that a particle of size A' is created, another particle of size A" = L{X') is also 
created. On the other hand, since the derivative of L is 

L'(A') = -(A-^ - A"i)(i-'*)/''A"'"\ (4.6) 
making the change of variables A" ^ X' = L^^(A"), we get 

Be{T,X,X")dX" ^ [ Se(^,A,(A'^-A"^)i/'^)(A'^-A"^)(i-'')/'^A"^"'dA'. (4.7) 
Hence, Be and Be are linked by 

Be{T, A, (A"^ - X"^f''^) {X^ - X'd){l'd)/dy,d-l ^ ^^^^^ 
or, since (A'' - }^d^(^-d)/dyd-l ^ _ y,<i)(d-l)/dy/i-d ^^^^ y, ^ ^^^^^ 

BeiT, A, (A"^ - X"dy/d^ (^d _ y/d)(l-d)/<iy/d-l ^ ^^^j^^ y ^4 g-j 

4.2 Building operator G 

Now, building an integral operator G to be used in (|2.5p . which takes those facts into account, 
consists in considering that the evolution of p(t, x. A) in the neighborhood of a given value of A, is 
the result of the following factors: a loss due to the aggregation of particles of size A with others, 
another loss due the fragmentation of particles of size A, a gain due to the aggregation of particles 
smaller than A and another gain due to the fragmentation of particles bigger than A. 

Quantifying the fragmentation-linked loss consists in noting that the density of particles at a 
given point x and in a given size A is nothing but p(t, x, A)/m(A), and in considering that, per unit 
of time, the number of particles of size A to fragment is in proportion with the number of present 
particles. Hence, 

^^BK^,A), (4.10) 
m(A) 

is the density, with respect to variables x and A, of particles of size A which fragment per unit of time 
in X and at t, when fluid conditions are J-. Then, the density of mass loss related to fragmentation 
is: 

"^(A) ^^g^ Bf{T, A) = Pit, X, A) BjiT, A). (4.11) 

In order to quantify loss linked to aggregation, it must be noted that the probability of a particle 
of size A to aggregate with a particle of size A', over a unit of time, is in direct proportion to the 
number of particles of size A and the number of particles of size A'. Consequently, the density, with 
respect to variables x and A, of particles of size A which aggregate is: 

^%^^%^^.(^,A,A')rfA', (4.12) 
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and the associated mass loss density, with respect to variables x and A, is 

p{t,^,X)B^h^Ba{T,\,X')dy. (4.13) 

The sum of the gain due to fragmentation is in two parts. The first part is the result of fragmen- 
tations, the smallest resulting particles of which are of size A, and the density of which, with respect 
to X and A can be written as: 

^(A) / ^^^^^ Bf{T, A') B,iT, A', A) dX'. (4.14) 

In order to understand (|4.14p . it has to be noted that p{t,x, X')/ni{X') is the density, with respect 
to variables x and A', of particles of size A'; Bf{J-^X') is the probability of a particle of size A' to 
fragment, within a unit of time; and, Be{J-,X' ,X) is the probability density function, with respect 
to variable A, of a particle of size A' to produce a particle of size A as its smallest resulting particle. 
Hence in (j4.14p . the integral is the density, with respect to variables x and A, of particles produced 
at size A as the smallest particule resulting from fragmentation. Multiplying this by m(A) gives the 
associated mass density. Because of (|4.4I) . (I4.14p may be written as 

^^^^p(M^ S/(^, A') B,{T, A', A) dX'. (4.15) 

A'eA ) 

The second part is the result of fragmentations whose largest resulting particles are of size A. 
The density, with respect to variables x and A, associated with this second part reads: 

^^^^ p(t,x AO ^^^^^ ^^^^^ ^^^g^ 
A'eA "^l^ ) 

or, because of (14. 9|) . 



A'eA m{X') 

As a consequence of (|4.15p and (|4.17p . it may be concluded that the mass-gain density, with 
respect to x and A, due to fragmentation result is 

A'eA m{>^) V / 

(4.18) 



The aggregation-related mass-gain density, with respect to x and A, reads 

'»<^)^1if^^^'^^ - A")"-"'A-X(^. A'. (A^ - A",./^, .A'. 

(4.19) 

or using (|4.2p . 

(4.20) 
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The form of the operator G which may be deduced from (|4.1ip . (|4.13p . (|4.18p and (|4.20p is 

G(^, p{t, X, .), A) = ~p{t, X, A) Bf{T, A) - / p{t, X, A) ^^^^^^ Ba{T, A, A') d\' 

(A)^i^li^ Bf{T, A') (Be{T, A', A) + BeiT, A', (A"^ - X'")'/'') (A"* - A'i)(i-'*)/''A'^-i) dX' 
m(X ) \ ) 



m 

A'eA 



(4.21) 

4.3 Properties of G 

If an operator G of the kind defined by (|4.2ip is chosen in equation p.Sp , then it does not influence 
the evolution of the total mass of sediment. In other words, a function p{t, A), not depending on x, 
solution to equation (|3.6p with G given by (|4.2ip . satisfies property (|3.7p since 

/ G(J^,p(t,.),A)dA = 0. (4.22) 

This may be seen by computing, on the one hand, the integral with respect to A of the third 
term of (|4.2ip . Using (|4.ip , it gives 

/ / M^)^^Bf{T,X')(B,iT,X',X) 

+ Be{T, A', (A"^ - X'^f'^) (A"^ - \d)i^-d)/dxd-T^^ dX' dX = 
I pit, A') BfiT, X')^Be{T, A', A) dA'dA 

p{t, X')Bf{T, X')Be{T, A', (A"^ - X'^Y''^) (A"^ - \d^(.^-d)/d^^d~i ^y^^. 
XeAJX'eA A 

(4.23) 

Making the change of variables (A, A') ^ (A, A')) with A = (A"* - X'^Y^'^ and A' = A', in the last 
integral, since X'^^^dX'dX = A'^^^dA'dA, it gives 

/ / p(t, X')BfiT, X')BeiT, A', A) A^^^^ A'^'^ dA'dA = 

J XeAJX'eA A 

f r yd _ \d 

/ / p(t,A')i3/(^,A')Se(^,A',A)— ^ dA'dA, (4.24) 
J XeAJX'eA 

and ((i?^ yields 

/ / MA)^^SK^,A')fs,(^,A',A) 
JxeAJx'eA n^(A') V 

+ Be{T, A', (A"* - A'^)!/'*) (A"' - A'i)(i-'i)/'iA''-i^ dA'dA = 

/ p{t,X')Bf{T,X')( [ BeiT,X',X) dx)dX' = [ p{t, X') Bf{T, X') dX' . (4.25) 
^a'gA \JxeA J Jx'eA 
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Once this computation is done, it is obvious that the integral with respect to A of the third term of 
(|4.2ip is the opposite of the integral with respect to A of the first term. 

On the other hand, integrating the last term of (I4.2ip with respect to A, gives 



2 



AeAJA'<A "i(A') m((A'' - A'")!/'!) 

(4.26) 



which, considering the following change of variables (A, A') i— > (A, A') = ((A"* — A''*)^/'^, A'), or equiv- 
alently (A, A') ^ (A, A') = ((A"* + A') for which X^-^dX'dX = X'^'^dX'dX, yields 



p{t,X)p{t,X') ' Ba{T,X\X)X'-''X''~'dX'dX = 

^ JxeAJx'eA miAjmiA ) 

\[ [ p{t,X)pit,X')(^ + ^—)BaiT,X',X)dX'dX = 
^JxeAJyeA \m{X) m(A')/ 

r r PM}pMlB^i:F,x',X)dx'dx, (4.27) 

AeAA'eA "^(A) 

where the first equality is obtained using (|4.ip . and the second using (|4.2p . As the last quantity in 
(|4.27p is nothing but the opposite of the integral, with respect to A, of the second term of (|4.2ip . 
(P:^ is true. 



5 On numerical applications 

In order to build numerical methods approximating equation (|2.5p . it must be noted that it is 
possible to make a splitting in time. For a given small time step At, this splitting routine consists, 
knowing an approximation p°(s, ., .) of p(s, ., .) at time s, in computing first p^{s + At, ., .), which is 
an approximation of the solution p to 



= 0, 



dx dy dz 

pis,;-) =P"(S,.,.), 

at time s + At, and then p°(s + At, ., .) as an approximation of the solution p to 

l^'^*^-"-^'- ,5.2) 

p{s,.,.)^ p\s + At,.,.), 

at time s + At. 

After discretizing A into / subsets A.;, an approximation of (|5.ip is made of a collection of / 
advection-diffusion equations (one for each A^), which are not mutually dependent. An approxi- 
mation of each of these advection-diffusion equations may be computed using any usual numerical 
advection-diffusion solver. Concerning the computation of an approximation of the solution to (|5.2p . 
which is the focus here, once the fluid field J- is known and the position in space is discretized, it 
comes down to computing a collection of approximations of p{s + At, A), solution to 

^^G{p,X), p{s,X)^po{X), (5.3) 
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(forgetting the dependence in T and x) for given functions po- 



In the case when G is given by (|4.2ip . one way to proceed would be to follow a Monte-Carlo 
method (see Lapeyre, Pardoux & Sentis [lOj) to approximate the integrals involved within the 
definition of G. Proceeding in this direction would also lead to numerical methods of the types used 
in Spilman & Levenspiel [15], Daniels Hughes [4], Meakin [12], Liffman [11], Shah et al. [14], Das 
[5], Spouge [16], Van Peborgh & Hounslow [171 and Kim et al. [9] 

There is another way in which an approximated solution of (|5.3|) with G given by (|4.2ip can be 
built. This way consists in building a discrete operator G from G, without breaking the structure 
yielding property (|4.22[) . As a matter of fact, a discrete version of property (|4.22p may be written 
for the discrete operator G. The method followed does have a relation to the method set out in 
Buet 2 , Rogier & Schneider [T3|, Degond & Lucquin 6J and Frenod & Lucquin 7 in the contexts 
of Boltzmann and Fokker-Plank equations. 

5.1 Discrete operator building 

First, as mentioned previously, A is discretized into / subsets A^ = [A^, A^+i) such that A^ n Aj = 
and A = uf^j^A^. Then, denoting by |Ai| the measure of A^, the following numbers are defined: 



B,}^—i S/(A)rfA, for^^l,...,/, (5.4) 
Iff I 



= WVTl / / — T77^a(A, A')dA'dA, for z and j = 1, . . . , /, (5.5) 

l^ill^jl J\eA^ Jx'eAj m(A') 



\Ai\\Aj\ JxeA, Jx'eA, to(A') 



Se(A', (A"* - A'*)i/'^)(A"^ - ~^d-^{i-d)/d~^d-i^~^,d _ ~^d-^ii-d)/d \ ^y^^^ f^^. ^ and j = 1, . . . ,/, (5.6) 



ga 



\^^\ JxeA. |Ajfc(A)| Jx'eA,,w m(A')m((A'^ - A"^)i/rf) 

BaiX, {X'^ ~ A"*)i/'^)(A'* - A"^)(i-'*'/''dA'dA, for ^ = 1, . . . , /, /c < i and ? < i, (5.7) 

where, in equation ((5?7)) . Ajfe(A) = ((A"* - Af^J^^'', (A'* - Af)^^''] n A^, which reads also Ajfe(A) = 
{A' e Afc, (A'^ - A"*)!/'' e A,}. The following functions are also defined: 

7 

B;/a) = 5]b,;ia,(a), 

i=l 

/ 7 

BJA, A') = E E B^:'' lA^W MA'), 

(5.8) 

b,/(a,a')-EEbJ1a.(a)1a,(a'), 

Bga(A, A',A*) = EEEBs'a' 1a,(A) lAj(A') l((A^-Af^i)i/d,(Ad-Af)i/d](A*), 
i=l k=\ l=\ 

where 1a stands for the characteristic function of set A.;. 
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With those definitions at hand, operator G{p{t, .), A) is set as 
G{p{t, .), A) = -p{t, A) By.(A) - / p{t, X)p{t, A') B,,(^, A, A') dX' 



A'eA 



+ / p{t,X')Bgf{X,X')dy + / p{t,X')p{t,X*)Bg^{X,X',X*)dX'dX*. (5.9) 

If p(t, .) is constant on every A,, with worth p^{t) or, in other words, if 

p(t,A) = ^p'(i) 1a.(A), (5.10) 

i=l 

then G{p{t, .), A) has the foUowing expression 

I II 

G(p(t,.),A) = -^B,y(t) 1a,(A) ~^^B,>Xi)P^W lA.(A) 

i=i i=i 'j=i 

II I i i 

+EEbJp''w iA.(A)+EEEB,f p'Wp'w ia.(a). (5.11) 

1=1 i=l i=l k=l 1=1 

5.2 Discrete operator properties 

By its construction, operator G{p{t, .), A) defined by (|5.9p is close to operator G{p{t, .), A) defined 
by l|4.2ip . as soon as p is regular enough. 

On the other hand, as a direct consequence of (|5.1ip . if p{t, .) is constant on every A.;, then 
G{p{t, .), A) is constant on every A^ and it is easy to see, as a consequence of the building of G, that 



G{p{t,.),X)dX^ / G(p(t,.),A)dA, (5.12) 
for i = 1, . . . , I. Then from equality (|4.22p . it may be deduced that 

G{p{t,.),X)dX^O. (5.13) 

AeA 

Hence it may be concluded that, if po is constant on every A^, the solution to 

^^G{p,X), p(s,A) = po(A), (5.14) 
is also constant on every A^ and satisfies 

/ p{t, A) dX is constant along time. (5.15) 

J\eA 

Consequently, a good way to build a mass-preserving numerical scheme to approximate ([5 
consists in approximating po by po defined by 



Po(A) = ^P^lA.(A) with pl^-^ f po{X)dX, (5.16) 
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and in approximating p{s + At, A), the solution to (|5.3p at time s + At, by p{s + At, A) defined as : 

p{.s + At, A) - po(A) + AtG{po{.), A). (5.17) 
The resuh p{s + At, A) wiU be close to p{s + At, A), constant on every A,; and will satisfy 

/ p{s + At,X)dX= [ po(A)dA = V|A,|p^= / po(A)dA. (5.18) 

6 Perspectives 

This paper puts forward a framework designed to process the characteristic evolution of sediment 
particles being transported by the water column. It gives simple examples of instantiations of this 
framework. 

It opens the way for many interesting questions. 

Concerning modeling, it would be useful to incorporate other sediment particle characteristics 
into the model, such as charge and fractal dimension, which seem important and are attentively 
studied by colloid scientists. To do this, new spaces A of larger dimension have to be built, taking 
into account the input from colloid sciences. 

Concerning mathematics, the existence of results for equations of the kind 

^^GiT,p,\), (6.1) 

with G given by l\'S.'6\i . (|3.1ip or (|4.2ip are an interesting challenge. 

Finally, concerning numerics, the path explored in section [5] has to be explored in greater depth, 
and software has to be designed to test the accuracy of such schemes. 

Acknowledgments - The author thanks Juliette Bouchery for proofreading the manuscript. 
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